﻿namespace GemmContest
{
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using MathNet.Numerics.Algorithms.LinearAlgebra;
    using MathNet.Numerics.LinearAlgebra.Double;
    using MathNet.Numerics.Threading;

    /// <summary>
    /// Used to time a GEMM implementation.
    /// </summary>
    /// <remarks>Submitters must use MathNet.Numerics.Threading.CommonParallel for any parallel code.</remarks>
    public class GemmBenchmark
    {
        /// <summary>
        /// Number of iterations
        /// </summary>
        private const int Iterations = 6;

        /// <summary>
        /// Matrix test sized.
        /// </summary>
        private static readonly int[] TestSizes = { 15, 32, 64, 99, 128, 129, 301, 512, 900, 1024, 1025 };

        /// <summary>
        /// The maximum number of threads to use for scalability testings.
        /// </summary>
        private static readonly int MaxThreads = Environment.ProcessorCount;

        /// <summary>
        /// Random number generator to fill the matrices - using a constant seed.
        /// </summary>
        private static readonly Random Random = new Random(0);

        /// <summary>
        /// Computes C := alpha * op(A) * op(B) + beta * C
        /// where op(X) is one of
        /// op(X) = X or op(X) = X'
        /// </summary>
        /// <param name="transposeA">To transpose matrix A or not.</param>
        /// <param name="transposeB">To transpose matrix B or not.</param>
        /// <param name="m">The number of rows of matrix op(A) and of the matrix C.</param>
        /// <param name="n">The number of columns of matrix op(B) and of the matrix C.</param>
        /// <param name="k">The number of columns of matrix op(A) and the number rows of the matrix op(B).</param>
        /// <param name="alpha">The value to scale the matrix A with.</param>
        /// <param name="a">The matrix A.</param>
        /// <param name="b">The matrix B.</param>
        /// <param name="beta">The value to scale the matrix B with.</param>
        /// <param name="c">The matrix C.</param>
        public static void Gemm(
            Transpose transposeA,
            Transpose transposeB,
            int m,
            int n,
            int k,
            double alpha,
            double[] a,
            double[] b,
            double beta,
            double[] c)
        {
            // Math.NET Numerics implementation
            // var rowsA = transposeA == Transpose.DontTranspose ? m : k;
            // var columnsA = transposeA == Transpose.DontTranspose ? k : m;
            // var rowsB = transposeB == Transpose.DontTranspose ? k : n;
            // var columnsB = transposeB == Transpose.DontTranspose ? n : k;
            // MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiplyWithUpdate(transposeA, transposeB, alpha, a, rowsA, columnsA, b, rowsB, columnsB, beta, c);
        }

        /// <summary>
        /// The main method.
        /// </summary>
        private static void Main()
        {
            // burn in
            RunSquareIteration(2);
            RunSquareIteration(2);

            var results = new List<Results>();

            for (var threadCount = 1; threadCount <= MaxThreads; threadCount++)
            {
                MathNet.Numerics.Control.NumberOfParallelWorkerThreads = threadCount;

                foreach (var order in TestSizes)
                {
                    Console.WriteLine("Testing: {0} thread(s) with {1} order matrices", threadCount, order);

                    var times = new List<double>();
                    for (var iteration = 0; iteration < Iterations; iteration++)
                    {
                        var time = RunSquareIteration(order); 
                        time += RunRectangularIteration(order);
                        times.Add(time);
                    }

                    // remove the fastest and slowest times
                    times.Sort();
                    times.RemoveAt(times.Count - 1);
                    times.RemoveAt(0);

                    results.Add(new Results { Threads = threadCount, Order = order, Time = times.Average() });
                }
            }

            var verySmall = BenchmarkVerySmall();
            var accuracy = ComputeAccuracy();
            DisplayResults(results, verySmall, accuracy);
        }

        /// <summary>
        /// Runs the iteration.
        /// </summary>
        /// <param name="order">The order.</param>
        /// <returns>The number of seconds the iteration took.</returns>
        private static double RunSquareIteration(int order)
        {
            var size = order * order;
            var a = new double[size];
            var b = new double[size];
            var c = new double[size];

            for (var index = 0; index < size; index++)
            {
                a[index] = Random.NextDouble();
                b[index] = Random.NextDouble();
                c[index] = Random.NextDouble();
            }

            var start = DateTime.Now;

            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, order, order, order, 0.0, a, b, 0.0, c);
            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, order, order, order, 1.0, a, b, 1.0, c);
            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, order, order, order, 0.0, a, b, 1.0, c);
            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, order, order, order, 1.0, a, b, 0.0, c);
            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, order, order, order, 2.2, a, b, 3.3, c);

            Gemm(Transpose.DontTranspose, Transpose.Transpose, order, order, order, 0.0, a, b, 0.0, c);
            Gemm(Transpose.DontTranspose, Transpose.Transpose, order, order, order, 1.0, a, b, 1.0, c);
            Gemm(Transpose.DontTranspose, Transpose.Transpose, order, order, order, 0.0, a, b, 1.0, c);
            Gemm(Transpose.DontTranspose, Transpose.Transpose, order, order, order, 1.0, a, b, 0.0, c);
            Gemm(Transpose.DontTranspose, Transpose.Transpose, order, order, order, 2.2, a, b, 3.3, c);

            Gemm(Transpose.Transpose, Transpose.DontTranspose, order, order, order, 0.0, a, b, 0.0, c);
            Gemm(Transpose.Transpose, Transpose.DontTranspose, order, order, order, 1.0, a, b, 1.0, c);
            Gemm(Transpose.Transpose, Transpose.DontTranspose, order, order, order, 0.0, a, b, 1.0, c);
            Gemm(Transpose.Transpose, Transpose.DontTranspose, order, order, order, 1.0, a, b, 0.0, c);
            Gemm(Transpose.Transpose, Transpose.DontTranspose, order, order, order, 2.2, a, b, 3.3, c);

            Gemm(Transpose.Transpose, Transpose.Transpose, order, order, order, 0.0, a, b, 0.0, c);
            Gemm(Transpose.Transpose, Transpose.Transpose, order, order, order, 1.0, a, b, 1.0, c);
            Gemm(Transpose.Transpose, Transpose.Transpose, order, order, order, 0.0, a, b, 1.0, c);
            Gemm(Transpose.Transpose, Transpose.Transpose, order, order, order, 1.0, a, b, 0.0, c);
            Gemm(Transpose.Transpose, Transpose.Transpose, order, order, order, 2.2, a, b, 3.3, c);

            return (DateTime.Now - start).TotalSeconds;
        }

        /// <summary>
        /// Runs the iteration.
        /// </summary>
        /// <param name="order">The order.</param>
        /// <returns>The number of seconds the iteration took.</returns>
        private static double RunRectangularIteration(int order)
        {
            const int Width = 4;
            var a = new double[order * Width];
            var b = new double[Width * order];
            var c = new double[order * order];

            for (var index = 0; index < a.Length; index++)
            {
                a[index] = Random.NextDouble();
            }

            for (var index = 0; index < b.Length; index++)
            {
                b[index] = Random.NextDouble();
            }

            for (var index = 0; index < c.Length; index++)
            {
                c[index] = Random.NextDouble();
            }

            var start = DateTime.Now;

            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, order, order, Width, 0.0, a, b, 0.0, c);
            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, order, order, Width, 1.0, a, b, 1.0, c);
            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, order, order, Width, 0.0, a, b, 1.0, c);
            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, order, order, Width, 1.0, a, b, 0.0, c);
            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, order, order, Width, 2.2, a, b, 3.3, c);

            Gemm(Transpose.DontTranspose, Transpose.Transpose, order, order, Width, 0.0, a, a, 0.0, c);
            Gemm(Transpose.DontTranspose, Transpose.Transpose, order, order, Width, 1.0, a, a, 1.0, c);
            Gemm(Transpose.DontTranspose, Transpose.Transpose, order, order, Width, 0.0, a, a, 1.0, c);
            Gemm(Transpose.DontTranspose, Transpose.Transpose, order, order, Width, 1.0, a, a, 0.0, c);
            Gemm(Transpose.DontTranspose, Transpose.Transpose, order, order, Width, 2.2, a, a, 3.3, c);

            Gemm(Transpose.Transpose, Transpose.DontTranspose, order, order, Width, 0.0, b, b, 0.0, c);
            Gemm(Transpose.Transpose, Transpose.DontTranspose, order, order, Width, 1.0, b, b, 1.0, c);
            Gemm(Transpose.Transpose, Transpose.DontTranspose, order, order, Width, 0.0, b, b, 1.0, c);
            Gemm(Transpose.Transpose, Transpose.DontTranspose, order, order, Width, 1.0, b, b, 0.0, c);
            Gemm(Transpose.Transpose, Transpose.DontTranspose, order, order, Width, 2.2, b, b, 3.3, c);

            Gemm(Transpose.Transpose, Transpose.Transpose, order, order, Width, 0.0, b, a, 0.0, c);
            Gemm(Transpose.Transpose, Transpose.Transpose, order, order, Width, 1.0, b, a, 1.0, c);
            Gemm(Transpose.Transpose, Transpose.Transpose, order, order, Width, 0.0, b, a, 1.0, c);
            Gemm(Transpose.Transpose, Transpose.Transpose, order, order, Width, 1.0, b, a, 0.0, c);
            Gemm(Transpose.Transpose, Transpose.Transpose, order, order, Width, 2.2, b, a, 3.3, c);

            return (DateTime.Now - start).TotalSeconds;
        }

        /// <summary>
        /// Benchmarks multiplying very small matrices.
        /// </summary>
        /// <returns>The time the benchmark took.</returns>
        private static double BenchmarkVerySmall()
        {
            const int Order = 3;
            const int Size = Order * Order;
            var a = new double[Size];
            var b = new double[Size];
            var c = new double[Size];

            for (var index = 0; index < Order; index++)
            {
                a[index] = Random.NextDouble();
                b[index] = Random.NextDouble();
                c[index] = Random.NextDouble();
            }

            const int Count = 250000;

            var start = DateTime.Now;
            for (var iteration = 0; iteration < Count; iteration++)
            {
                Gemm(Transpose.DontTranspose, Transpose.DontTranspose, Order, Order, Order, 0.0, a, b, 0.0, c);
                Gemm(Transpose.DontTranspose, Transpose.DontTranspose, Order, Order, Order, 1.0, a, b, 1.0, c);
                Gemm(Transpose.DontTranspose, Transpose.DontTranspose, Order, Order, Order, 0.0, a, b, 1.0, c);
                Gemm(Transpose.DontTranspose, Transpose.DontTranspose, Order, Order, Order, 1.0, a, b, 0.0, c);
                Gemm(Transpose.DontTranspose, Transpose.DontTranspose, Order, Order, Order, 2.2, a, b, 3.3, c);

                Gemm(Transpose.DontTranspose, Transpose.Transpose, Order, Order, Order, 0.0, a, b, 0.0, c);
                Gemm(Transpose.DontTranspose, Transpose.Transpose, Order, Order, Order, 1.0, a, b, 1.0, c);
                Gemm(Transpose.DontTranspose, Transpose.Transpose, Order, Order, Order, 0.0, a, b, 1.0, c);
                Gemm(Transpose.DontTranspose, Transpose.Transpose, Order, Order, Order, 1.0, a, b, 0.0, c);
                Gemm(Transpose.DontTranspose, Transpose.Transpose, Order, Order, Order, 2.2, a, b, 3.3, c);

                Gemm(Transpose.Transpose, Transpose.DontTranspose, Order, Order, Order, 0.0, a, b, 0.0, c);
                Gemm(Transpose.Transpose, Transpose.DontTranspose, Order, Order, Order, 1.0, a, b, 1.0, c);
                Gemm(Transpose.Transpose, Transpose.DontTranspose, Order, Order, Order, 0.0, a, b, 1.0, c);
                Gemm(Transpose.Transpose, Transpose.DontTranspose, Order, Order, Order, 1.0, a, b, 0.0, c);
                Gemm(Transpose.Transpose, Transpose.DontTranspose, Order, Order, Order, 2.2, a, b, 3.3, c);

                Gemm(Transpose.Transpose, Transpose.Transpose, Order, Order, Order, 0.0, a, b, 0.0, c);
                Gemm(Transpose.Transpose, Transpose.Transpose, Order, Order, Order, 1.0, a, b, 1.0, c);
                Gemm(Transpose.Transpose, Transpose.Transpose, Order, Order, Order, 0.0, a, b, 1.0, c);
                Gemm(Transpose.Transpose, Transpose.Transpose, Order, Order, Order, 1.0, a, b, 0.0, c);
                Gemm(Transpose.Transpose, Transpose.Transpose, Order, Order, Order, 2.2, a, b, 3.3, c);
            }

            return (DateTime.Now - start).TotalSeconds;
        }

        /// <summary>
        /// Displays the results of the benchmark.
        /// </summary>
        /// <param name="results">Contains the result of each run.</param>
        /// <param name="verySmall">The time of the very small benchmark.</param>
        /// <param name="accuracy">Accuracy values of each run.</param>
        private static void DisplayResults(IEnumerable<Results> results, double verySmall, Tuple<double, double> accuracy)
        {
            var total = results.Sum(item => item.Time);

            Console.WriteLine("Threads, Order, Average Time");
            foreach (var result in results)
            {
                Console.WriteLine("{0}, {1}, {2}", result.Threads, result.Order, result.Time);
            }

            Console.WriteLine();
            Console.WriteLine("Sum of all averaged times: {0}", total);
            Console.WriteLine("Very small time: {0}", verySmall);
            Console.WriteLine("Maximum Absolute Error: {0}", accuracy.Item1);
            Console.WriteLine("Root Mean Square Error: {0}", accuracy.Item2);
        }

        /// <summary>
        /// Computes the accuracy of the GEMM implementation.
        /// </summary>
        /// <returns>The maximum absolute error and the root mean square error.</returns>
        private static Tuple<double, double> ComputeAccuracy()
        {
            var matlabReader = new MathNet.Numerics.LinearAlgebra.IO.MatlabMatrixReader<double>("gemm.mat");
            var a = (DenseMatrix)matlabReader.ReadMatrix("a");
            var b = (DenseMatrix)matlabReader.ReadMatrix("b");
            var c = (DenseMatrix)matlabReader.ReadMatrix("c");
            var testValues1 = (DenseMatrix)matlabReader.ReadMatrix("d");
            var testValues2 = (DenseMatrix)matlabReader.ReadMatrix("e");
            var testValues3 = (DenseMatrix)matlabReader.ReadMatrix("f");
            var testValues4 = (DenseMatrix)matlabReader.ReadMatrix("g");

            var max = 0.0;
            var rmse = 0.0;
            var count = c.RowCount * c.ColumnCount;

            var testC = (DenseMatrix)c.Clone();
            Gemm(Transpose.DontTranspose, Transpose.DontTranspose, a.RowCount, b.ColumnCount, a.ColumnCount, 2.0, a.Data, b.Data, 3.0, testC.Data);
            for (var index = 0; index < count; index++)
            {
                var diff = testValues1.Data[index] - testC.Data[index];
                max = Math.Max(max, Math.Abs(diff));
                rmse = diff * diff;
            }

            testC = (DenseMatrix)c.Clone();
            Gemm(Transpose.DontTranspose, Transpose.Transpose, a.RowCount, b.ColumnCount, a.ColumnCount, 2.0, a.Data, a.Data, 3.0, testC.Data);
            for (var index = 0; index < count; index++)
            {
                var diff = testValues2.Data[index] - testC.Data[index];
                max = Math.Max(max, Math.Abs(diff));
                rmse = diff * diff;
            }

            testC = (DenseMatrix)c.Clone();
            Gemm(Transpose.Transpose, Transpose.DontTranspose, a.RowCount, b.ColumnCount, a.ColumnCount, 2.0, b.Data, b.Data, 3.0, testC.Data);
            for (var index = 0; index < count; index++)
            {
                var diff = testValues3.Data[index] - testC.Data[index];
                max = Math.Max(max, Math.Abs(diff));
                rmse = diff * diff;
            }

            testC = (DenseMatrix)c.Clone();
            Gemm(Transpose.Transpose, Transpose.Transpose, a.RowCount, b.ColumnCount, a.ColumnCount, 2.0, b.Data, a.Data, 3.0, testC.Data);
            for (var index = 0; index < count; index++)
            {
                var diff = testValues4.Data[index] - testC.Data[index];
                max = Math.Max(max, Math.Abs(diff));
                rmse = diff * diff;
            }

            rmse = Math.Sqrt(rmse / (count * 4));

            return new Tuple<double, double>(max, rmse);
        }

        /// <summary>
        /// Object to hold timing results.
        /// </summary>
        private struct Results
        {
            /// <summary>
            /// The number of threads used in a run.
            /// </summary>
            public int Threads;

            /// <summary>
            /// The size of the matrix in a run.
            /// </summary>
            public int Order;

            /// <summary>
            /// The average time it took to complete a run.
            /// </summary>
            public double Time;
        }
    }
}